//
// GPSMapEdit
// (c) Konstantin Galichsky (kg@geopainting.com), 2002-2004
//
// New Zealand Map Grid projection.
//

# include "StdAfx.h"
# include "NewZealandMapGrid.h"
# include "PtAlgo.h"
# include "Globals.h"

//# include <complex>
class complex_t {
  public:
	double re, im;

	complex_t (double _re, double _im) : re (_re), im (_im) {}
	complex_t (const complex_t & _other) : re (_other.re), im (_other.im) {}

	double real () const {return re;}
	double imag () const {return im;}

	void operator += (const complex_t & _other) {
		re += _other.re;
		im += _other.im;
	}
	void operator *= (const complex_t & _other) {
		const double r = re*_other.re - im*_other.im;
		const double i = im*_other.re + re*_other.im;
		re = r;
		im = i;
	}
	void operator *= (double _other) {
		re *= _other;
		im *= _other;
	}
};

complex_t operator * (const complex_t & _1, int _2) {
	complex_t c (_1);
	c *= _2;
	return c;
}

complex_t operator / (const complex_t & _1, const complex_t & _2) {
	const double a2 = _2.re*_2.re + _2.im*_2.im;
	complex_t c (
		(_1.re*_2.re + _1.im*_2.im)/a2,
		(_1.im*_2.re - _1.re*_2.im)/a2
	);
	return c;
}

//////////////////////////////////////////////////

	// International ellipsoid (Geodetic Datum 1949)
const double c_fNZMG_a = 6378388.0;      // Semi-major axis, in meters

	// NZMG projection parameters
const double c_fOriginLat = -41.0;       // Latitude of origin, in degrees
const double c_fOriginLon = 173.0;       // Longitude of origin, in degrees
const double c_fFalseNorthing = 6023150; // False northing, in meters
const double c_fFalseEasting  = 2510000; // False easting, in meters

	// Minimum/maximum easting/northing values
const double c_fMinEasting  = 1810000;
const double c_fMaxEasting  = 3170000;
const double c_fMinNorthing = 5160000;
const double c_fMaxNorthing = 6900000;

const double c_fMaxLat = -33.5;
const double c_fMinLat = -48.5;
const double c_fMaxLon = 180.0;
const double c_fMinLon = 165.5;

const complex_t c_B [6] = {
	complex_t (0.7557853228, 0.0 ),
	complex_t (0.249204646,  0.003371507),
	complex_t (-0.001541739, 0.041058560),
	complex_t (-0.10162907,  0.01727609 ),
	complex_t (-0.26623489,  -0.36249218),
	complex_t (-0.6870983,   -1.1651967 )
};

const complex_t c_C [6] = {
	complex_t (1.3231270439, 0.0         ),
	complex_t (-0.577245789, -0.007809598),
	complex_t (0.508307513,  -0.112208952),
	complex_t (-0.15094762,  0.18200602  ),
	complex_t (1.01418179,   1.64497696  ),
	complex_t (1.9660549,    2.5127645   )
};

const double c_D [9] = {
	1.5627014243, 0.5185406398, -0.03333098, -0.1052906, -0.0368594,
	0.007317,     0.01220,      0.00394,     -0.0013
};

void CNewZealandMapGrid::ApplyTo (point_t & _point) const {
//	assert (c_fMinEasting  <= _point.x && _point.x <= c_fMaxEasting);
//	assert (c_fMinNorthing <= _point.y && _point.y <= c_fMaxNorthing);

	const complex_t z (
		(_point.y - c_fFalseNorthing)/c_fNZMG_a,
		(_point.x - c_fFalseEasting )/c_fNZMG_a
	);

	complex_t Zeta (c_C [5].real (), c_C [5].imag ());
	int n;
	for (n = 4; n >= 0; -- n) {
		Zeta *= z;
		Zeta += c_C [n];
	}
	Zeta *= z;

	for (size_t i = 0; i < 2; ++ i) {
		complex_t Numer = c_B [5];
		Numer *= 5.;
		complex_t Denom = c_B [5];
		Denom *= 6.;

		for (n = 4; n >= 1; -- n) {
			Numer *= Zeta;
			Numer += c_B [n]*n;

			Denom *= Zeta;
			Denom += c_B [n]*(n + 1);
		}

		Numer *= Zeta;
		Numer *= Zeta;
		Numer += z;

		Denom *= Zeta;
		Denom += c_B [0];

		Zeta = Numer/Denom;
	}

	double dphi = 0;
	for (n = 8; n >= 0; -- n) {
		dphi += c_D [n];
		dphi *= Zeta.real ();
	}
	_point.y = c_fOriginLat + (dphi*1e5/3600);
//	assert (c_fMinLat <= _point.y && _point.y <= c_fMaxLat);

	_point.x = c_fOriginLon + Zeta.imag ()*180/c_PI;
//	assert (c_fMinLon <= _point.x && _point.x <= c_fMaxLon);
}

///////////////////////////////

# ifdef _DEBUG
void Test_NZMG () {
	point_t p (2693473.80f, 6687001.06f);
	CNewZealandMapGrid nzmg;
	nzmg.ApplyTo (p);
	assert (::fabs (p.x - 175.0f) < 1e-6);
	assert (::fabs (p.y - -35.0f) < 1e-6);
}
# endif // _DEBUG
